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ABSTRACT 

Axisymmetric magnetorotational instability (MRI) in viscous accretion disks is investigated by 
linear analysis and two-dimensional nonlinear simulations. The linear growth of the viscous MRI is 
characterized by the Reynolds number defined as -Rmri = where va is the Alfven velocity, 

V is the kinematic viscosity, and is the angular velocity of the disk. Although the linear growth 
rate is suppressed considerably as the Reynolds number decreases, the nonlinear behavior is found to 
be almost independent of i?MRi- At the nonlinear evolutionary stage, a two-channel flow continues 
growing and the Maxwell stress increases until the end of calculations even though the Reynolds 
number is much smaller than unity. A large portion of the injected energy to the system is converted 
to the magnetic energy. The gain rate of the thermal energy, on the other hand, is found to be 
much larger than the viscous heating rate. Nonlinear behavior of the MRI in the viscous regime 
and its difference from that in the highly resistive regime can be explained schematically by using the 
characteristics of the linear dispersion relation. Applying our results to the case with both the viscosity 
and resistivity, it is anticipated that the critical value of the Lundquist number S'mri = v\/rjVl for 
active turbulence depends on the magnetic Prandtl number S'mri.c oc Pm}^^ in the regime of Pm ^ 1 
and remains constant when Pm <^ 1, where Pm = S'iviriZ-Rmri = ^ Iv ^-nd rj is the magnetic diffusivity. 
Subject headings: accretion, accretion disks — MHD — turbulence — methods: numerical 



1. INTRODUCTION 

Magnetohydrodynamic (MHD) turbulence is the most 
promising candidate for angular momentum transport in 
astrophysical disk systems. Nonlinear behaviors of the 
magnetorotational instability (MRI) are actively investi- 
gated as a driving mechanism of MHD turbulence over 
the last decades (Balbus & Hawley 1991, 1998). The 
central issue in MRI research is its nonlinear properties, 
in particular, the saturation amplitude of the instabil- 
ity. Although the key processes governing the nonlin- 
ear saturation are scoped by global and local numerical 
studies, it is not fully explained yet (Hawley & Balbus 
1992; Hawley et al. 1995, 1996; Brandenburg et al. 1995; 
Matsumoto & Tajima 1995; Stone et al. 1996; Hawley 
' 2000; Machida et al. 2000; Arlt & Riidiger 2001; Bal- 
; bus 2003). Recently, Lesur & Ogilvie (2008) argue that, 
in the shearing box context with zero-net vertical flux, 
MHD turbulence is sustained through nonlinear classical 
dynamo activity once the MRI is operated. It would 
be necessary and significative to study the saturation 
process from the microscopic viewpoint of the physical 
sustaining mechanism for MHD turbulence as they have 
done. 

Ohmic dissipation is one of the crucial processes that 
determine the saturation amplitude of the MRI. Linear 
growth rate of the MRI can be reduced significantly be- 
cause of the suppression by ohmic dissipation. Two- and 
three-dimensional local shearing box simulations (Sano 
et al. 1998, 2004; Sano & Stone 2002) have shown that 
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physical properties of the saturated turbulence depend 
on the Lundquist number S'mri = v\/riQ., where va is 
the Alfven velocity, is the angular velocity, and -q is 
the magnetic diffusivity (see also Fleming et al. 2000; 
Ziegler & Riidiger 2001; Liu et al. 2006). Particularly, 
when 5'mri *C 1, the size of the saturated stress decreases 
with the decrease of S'mri (Sano & Stone 2002; Pessah 
et al. 2007). It is also pointed out that magnetic recon- 
nection plays an important role in the energy dissipation 
of MRI driven turbulence (Sano & Inutsuka 2001). 

Numerical issues are one of the main reasons why the 
saturation physics of the MRI is remained to be under- 
stood. Fromang & Papaloizou (2007) demonstrate the 
efficiency of angular momentum transport decreases lin- 
early with the grid spacing as the resolution increases. 
Although it is very difficult to distinguish between the 
numerical and physical factors working as the saturation 
mechanism (King et al. 2007; Silvers 2007), current re- 
searches of the MRI pay much attention to the numerical 
factors with the greatest care. Pessah et al. (2007) de- 
rive a scaling law of the saturated stress from past wide 
variety of numerical results by analytically eliminating 
the numerical factors such as box size and resolutions. 

More straightforward way for decontaminating the nu- 
merical factors is to bring explicitly the physical diffusiv- 
ities much larger than the numerical one into the com- 
putational study. Lesur & Longaretti (2007) have per- 
formed first systematic study of the MRI in the pres- 
ence of both viscous and magnetic dissipations, which 
are larger than the numerical diffusivities. For the cases 
with nonzero net flux of the vertical field, the transport 
property in the saturated state depends on the magnetic 
Prandtl number Pm = v/r], where v is the kinematic 
viscosity. Fromang et al. (2007) have reported that in 
zero magnetic flux cases the turbulent activity is an in- 
creasing function of the magnetic Prandtl number Pm. 
Linear behaviors of the MRI in the presence of both the 
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viscosity and resistivity are analytically studied by Pes- 
sah & Chan (2008) comprehensively. 

The magnetic Prandtl number Pm takes a wide range 
of values in astrophysical disk systems. In protoplan- 
etary disks surrounding young stellar objects, the mag- 
netic Prandtl number is much smaller than unity because 
of their low ionization degree (Nakano 1984; Umebayashi 
& Nakano 1988; Sano et al. 2000). In accretion disks 
of compact X-ray sources and active galactic nuclei, the 
magnetic Prandtl number ranges from ~ 10"'^ to 10'^ 
depending on the distance from the central object (Bal- 
bus & Henri 2008). In coUapsar disks which is known as 
the central engine of gamma-ray bursts (Woosley 1993), 
the physical state with high magnetic Prandtl number 
of Pm > 10^" is expected to be realized in their evolu- 
tionary stage as a result of the large neutrino viscosity 
(Masada et al. 2007). Therefore, more systematic and 
deeper study on the MRI in the presence of both the vis- 
cosity and resistivity is quite important for understand- 
ing the accretion process triggered by the MRI in various 
disk systems. 

One important unsettled matter, in these situations, 
is the role of the kinematic viscosity at the nonlinear 
stage of the MRI. In general, the viscosity as well as 
the magnetic resistivity can suppress the growth of the 
MRI. However the dependence of nonlinear outcome on 
the Prandtl number indicates that the role of the vis- 
cosity in MRI turbulence could be different from that 
of the resistivity. Then, the main purpose of this pa- 
per is to reveal nonlinear features of the MRI in vis- 
cous accretion disks. As is described in what follows, lin- 
ear growth of the MRI is characterized by the Reynolds 
number i?MRi = v^/vfl in the viscous fluid, and by the 
Lundquist number Sum = v^/rifl in the resistive fluid. 
Focusing on these two non-dimensional parameters, we 
clarify the difference in nonlinear behaviors of the MRI 
between the viscous and resistive systems. 

Our paper is organized as follows. In § 2, linear fea- 
tures of the MRI in the viscous fluid are presented. In 
§ 3, nonlinear behavior of the MRI is investigated by 
two-dimensional MHD simulations taking account of the 
viscous terms. The differences between the effect of the 
viscosity and resistivity are also clarified in § 3. Finally 
we make an physical explanation for our nonlinear re- 
sults with the help of the linear dispersion relation. Ap- 
plying our results to double diffusive systems, we predict 
a condition for sustaining active MRI turbulence in the 
presence of both the viscosity and resistivity in § 4. 

2. LINEAR ANALYSIS 

First, we provide the linear features of the MRI in a 
viscous accretion disk threaded by a uniform vertical field 
Bz- Plane- wave perturbation theory, with WKB spatial 
and temporal dependence cx eyip{ikzZ + "ft), gives a local 
axisymmetric dispersion equation for the MRI, 
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(Menou et al. 2006; Masada et al. 2007; Pessah & Chan 
2008), where 7 = 7/SI is the growth rate normalized 
by angular velocity fJ, and k^ — k^VA/^ is the verti- 



cal wavenumber normalized by fl/vA- The epicyclic fre- 
quency, normalized by the angular velocity, can be ex- 
pressed as K = [2(2 — q)]^/^ by using the shear parameter 
q = — dlnfi/dlnr. The Reynolds number for the MRI 
is defined as i?MRi — VL/u = v\/v^. Here the charac- 
teristic velocity and length are V = va and L = va/^, 
respectively. In this paper, we focus on the MRI in Kep- 
lerian disks where the epicyclic frequency is equal to the 
angular velocity {q = 3/2). 

The dispersion equation ([T]) is characterized by the 
Reynolds number i?MRi. The linear growth rate of the 
MRI is shown as a function of the vertical wavenumber 
for the cases i?MRi = 0.1, 1.0, 10.0 and 00 in Figure [TJ 
When the kinematic viscosity is negligible (i?MRi ^ 1), 
the dispersion relation is identical to that of the ideal 
MHD case. If the Reynolds number is less than unity, 
the growth of the MRI is suppressed and the maximum 
growth rate is reduced significantly. The most unstable 
wavenumber decreases with decreasing i?MRi, because 
the viscous damping becomes more efficient for shorter 
wavelength perturbations. However, it is interesting that 
the critical wavenumber for the instability, fc^^crit — \/2q, 
remains unchanged in spite of the size of -Rmri- 

Figures [2^ and[2lD show the growth rate and wavenum- 
ber of the fastest growing mode as a function of the 

Reynolds number. These figures indicate that 7niax and 
~ 1/2 

fc2,max are proportional to Ryi^i when i?MRi ^ 1 (Pes- 
sah & Chan 2008). In the regime of i?MRi ^ 1, the 
dispersion equation (U)) can be simplified and reduced to 



[K^ + 2(2 - q)] - 2qK'^ = , 
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where G = 7/-Rmri and K = The fastest grow- 

ing wavelength and the maximum growth rate are ob- 
tained analytically from this equation; 
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2(2-9). 

Qualitative features of the fastest growing mode can 
be explained as follows: Keplerian shear flow is a key 
ingredient of the unstable growth of the MRI. The vis- 
cous dissipation affects the growth of the MRI when the 
damping rate is comparable to the shear rate for a per- 
turbation. Then the viscous damping rate of the unstable 
mode would be balanced with the shear rate in the vis- 
cous regime, k'^v ^ |df2/dlnr| ~ VI. This relation gives 

the most unstable wavenumber kz^max — {^/vY^'^, and 

1/2 ' 
thus fcz^inax oc Ryiiii- Siuce the maximum growth rate of 

the MRI is equal to the Alfven frequency of the fastest 
growing mode, it corresponds to 7max — fcz.max'^A) or 
1/2 

7max OC RyiYii- This is why the normalized wavenumber 
and growth rate of the fastest growing mode are pro- 
portional to the i?MRi presence of large viscous 
dissipation. 

Note that the i?MRi-dependence of the fastest grow- 
ing mode is slightly different from that on the Lundquist 
number (S'mri = 'l>a/v^) in presence of ohniic dissi- 
pation (Sano & Miyama 1999), where 77 is the magnetic 
diffusivity. Based on the linear analysis, the resistiv- 
ity can suppress the MRI more efficiently compared to 
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the viscosity. In the resistive regime, the ohmic dissi- 
pation can suppress the MRI when the dissipation time 
of the unstable mode is comparable to the Alfven time, 
^rnax/v — ^max/vA, Thus the most Unstable wavelength 
is given by A^ax — fl/'^^A- The growth rate of the fastest 
growing mode is given by 7inax ^ ■^A/Amax ^ Va/v- 
Then we can obtain the relations /c^.max oc S'mri and 
7max oc jS'mri- It is Stressed that the critical wavenumber 
fcz.crit in the resistive regime is also proportional to the 
Lundquist number, k^.^a oc S'mri- Unstable wavelength 
and growth rate of the MRI for ideal MHD, resistive, and 
viscous cases are summarized in Table [T] 

3. NONLINEAR ANALYSIS 

3.1. Numerical Setting 

For elucidating the nonlinear features of the MRI, vis- 
cous MHD equations are solved with a finite-differencing 
code which was developed by Sano et al. (1998). The hy- 
drodynamic module of our scheme is based on the second- 
order Godunov scheme (van Leer 1979), which consists 
of Lagrangian and remap steps. The Riemann solver is 
modified for accounting the effect of tangential magnetic 
fields. The field evolution is calculated with the Consis- 
tent MoC-CT method (Clarke 1996). The energy equa- 
tion is solved in the conservative form and the viscous 
terms are calculated in the Lagrangian step. The advan- 
tages of our scheme are its robustness for strong shocks 
and the satisfaction of the divergence-free constraint of 
magnetic fields (Evans & Hawley 1988; Stone & Norman 
1992). 

We use a local shearing box model (Hawley et al. 
1995). In this approximation, equations of viscous MHD 
are written in a local Cartesian frame of reference coro- 
tating with the disk at the angular velocity f2 corre- 
sponding to a fiducial radius R. Then the coordinates 
are presented as x = r — R, y = R<j) — VLt, and z. The 
fundamental equations are written in terms of these co- 
ordinates within a small region surrounding the fiducial 
radius, in Ar <C i?. 
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where e is the specific internal energy, ^ 
and X is the bulk viscosity. In the following, the effect 
of the bulk viscosity is neglected (x <C i^)- The term 



2qfl^x in the momentum equation is the tidal expansion 
of the effective potential. Assuming the ideal gas, the 
pressure is given by P = (7— l)pe. The spatially uniform 
kinematic viscosity and a constant ratio of specific heats 
(7 — 5/3) are considered. 

Since we focus on the local properties of the instability, 
we employ a numerical grid representing a small section 
of the disk interior for a local disk model. Adopting the 
Keplerian shear flow as unperturbed state, the azimuthal 
velocity is given by Vy = —qQx in the frame corotating 
with the velocity R^l. The initial field geometry is a weak 
uniform field in the vertical direction Bz = Bq. The 
radial force balance at the initial state is thus realized 
between the Coriolis force and the tidal force. 

Two-dimensional calculation is performed in the 
radial-vertical plane with a volume bounded by a; = z = 
±H/2, where H = {2/^Y/'^Cs/Vl is the scale height of 
the disk. We use a uniform grid of 128 x 128 zones. A 
periodic boundary condition is applied in the vertical di- 
rection. For the radial boundary condition, we adopt 
a sheared periodic boundary condition (Hawley et al. 
1995). In this model, the vertical component of grav- 
ity can be ignored. Except for the shear velocity, the 
physical quantities are thus assumed to be spatially uni- 
form; p = pq and P = Pq where po and Pq are con- 
stant values. We choose normalizations with po = 1; 

= 1, f7 = IQ-^ and Pq = 5 x 10"^ Initial pertur- 
bations are introduced as spatially uncorrelated velocity 
and adiabatic pressure fluctuations. These fluctuations 
have a zero mean value with a maximum amplitude of 
\6P\/Pq = 10-2 and \dv\/cs = lO'^. 

Our local disk model is characterized by non- 
dimensional parameters, which are the plasma beta of 
the initial field strength (3o = SttPq/Bq and the initial 
Reynolds number Pmrj. In this paper, we show the re- 
sults focusing on the effects of the Reynolds number. In 
what follows, we fix the initial field strength as /3o = 10* 
in all our models. In the case Pmri ^ 1, the fastest 
growing wavelength of the local disk system Avis can be 
defined from equation ([5]), 

Amax 27r/_t/\l/2 2tT 



Av 



H H \ni vm^i' ^^^^ 

where Amax is the fastest growing wavelength expected 
by the linear analysis for viscous MHD case. To capture 
the most unstable mode in the computational domain, 
-^vis S 1> the Reynolds number must be Pmri ^ 0.004. 

3.2. Results 
3.2.1. Growth Rate at ttie Linear Phase 

We investigate the growth rate at the early linear phase 
in order to confirm the accuracy of our numerical scheme. 
All the simulations begin with random perturbation of 
very small amplitude so that any growing modes should 
be well described by a linear analysis during the first few 
orbits of the evolution. The time history of each mode 
is followed through a two-dimensional Fourier decompo- 
sition carried out at frequent time intervals. We define 
the Fourier coefficient ak of the radial velocity Vx as 

«fc(^) = 77T-T2 / i Vx{x,z,t)e'^'''dxdz . (13) 



(27r)2 J , 

Here we select the modes with h 
the most unstable. 



because they are 
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The numerical growth rates in early linear phase are 
plotted over the analytic dispersion relation in Figure [TJ 
The growth rate at = 0.32, 0.64, 0.97, and 1.28 are 
shown as representative cases, which are obviously re- 
producing the analytic results for all the cases. Thus our 
numerical scheme has ability to simulate correctly the 
unstable growth of perturbations excited by the MRI in 
the presence of the viscous dissipation. 

3.2.2. Dependence on the Reynolds Number 

The efficiency of angular momentum transport is given 
by the x-y component of the stress tensor. 



Wxy = WM + Wr 



BxBy 

An 



+ pVxSvy , (14) 



where wm and wr are Maxwell and Reynolds stresses, re- 
spectively. This is related to the a parameter of Shakura 
& Sunyaev (1973) by a = w^y/P = {wm + wr)/P. The 
Maxwell stress is proportional to the magnetic energy 
and usually dominates over the Reynolds stress in MRI 
driven turbulence. 

To demonstrate the nonlinear features of the MRI, 
the time-evolutions of volume-averaged Reynolds and 
Maxwell stresses, {ur) = {wr)/{P) and (aA/) = 
{w]\i)/{P), are depicted in Figure [3] for the cases with 
different Reynolds numbers i?MRi = 0.01, 0.1, 1.0, and 
10.0. The single bracket indicates a volume average of 
physical quantities. The horizontal axis is the time nor- 
malized by the rotation time trot = 27r/f2. The kinematic 
viscosity in those models is equivalent to the a parameter 
of the size ~ {RmriPo)~^ ^ 0.01. 

The linear growth rate decreases as the Reynolds num- 
ber decreases. After the linear growth of the MRI, a two- 
channel flow appears for the ideal MHD cases (Hawley 
& Balbus 1992). The two-channel flow is an axisymmet- 
ric MRI mode whose vertical wavelength fits the vertical 
box size. This linearly unstable mode is also an exact so- 
lution of nonlinear MHD equations, so that the magnetic 
field can be amplified exponentially even at the nonlin- 
ear regime (Goodman & Xu 1994). Similar behavior is 
found in all the viscous models even though the Reynolds 
number is much smaller than unity. The magnetic energy 
continues growing and is not saturated even at the non- 
linear regime. The Maxwell stress increases until the end 
of calculations for all the models. The Reynolds stress, 
on the other hand, approaches a constant value at the 
nonlinear stage. 

The time evolution of (atot), which is the sum of {(xr) 
and (om), is shown in Figure |4^. The nonlinear behavior 
of the MRI in viscous fluid is quite different from that 
in the models taking account of the magnetic diffusivity 
(Sano et al. 1998, 2004; Fleming et al. 2000; Sano & In- 
utsuka 2001). For the purpose of comparison, Figure[4jD 
shows (tttot) in the resistive MHD runs with different 
Lundquist numbers S'mri- The initial conditions are the 
same as the viscous models except for the dissipation 
terms. When the dissipation processes can be negligible 
(-Rmri ^ 1 and S'mri ^ 1), the evolution is quite sim- 
ilar to that of ideal MHD. The viscous dissipation can- 
not damp MRI driven turbulence in the nonlinear regime 
even though i?MRi <C 1. In contrast, the MRI saturates 
and MHD turbulence dies away at the nonlinear regime 
for the case with 5mri ^ 1. 



3.2.3. Energy Injection mto the Local System 
The total energy within the shearing box is defined as 



T = JdV[p{v^/2 + € + ^)+B^ 



where 



is the tidal expansion of the effective potential (Hawley 
et al. 1995). Using the evolution equations for viscous 
MHD system [eqs. ([5|)-([TT|)]. the time derivative of the 
total energy gives 
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where dA is the surface element and the integral is taken 
over either of the radial boundaries. We can derive the 
last term of above equation by assuming that the radial 
gradient of the perturbed azimuthal velocity is negligi- 
ble, that is dvy/dx ~ — gfJ. The energy injection rate 
through the radial boundary is proportional to the kine- 
matic viscous stress as well as the turbulent stress Wxy 
at the boundary. 

Using the volume-averaged values instead of the 
surface-averaged ones at the radial boundary, the vol- 
ume average of the total energy changing rate (i?tot) is 
given by 

(.Btot) = {hn.,^) + (^i„,„) , (16) 

where {Ein,w) = q_^{wxy) is the energy injection rate 
caused by the turbulent stress (Sano & Inutsuka 2001), 
and {Ein,v) = vq^^^{p) is done by the kinematic viscous 
stress. The second term is always positive. When the 
turbulent stress is positive, the total energy of the sys- 
tem must increase. The source of the input energy is 
the background shear motion. In realistic disk systems, 
positive stresses lead to inward mass accretion, bringing 
a loss of gravitational energy. The gain in total energy 
in the shearing box represents this energy release. 

The energy budget in our simulations can be satisfying 
the relation (|16p. because our numerical scheme solves 
the energy equation in terms of the total energy. Fig- 
ure [5] shows the time evolution of (Stot), {Em^w) and 
{Ein.v) for the case with i?MRi = 0.01. The vertical axis 
is given in the unit of -Etho/^rot, where ii^tho = Po/ij-l) 
is the initial thermal energy. During the early linear 
phase until about 15 orbits, the kinematic viscous stress 
takes a major role in the energy injection. In contrast, 
the turbulent stress becomes predominant after the MRI 
grows sufficiently. The sum of these contributions is ex- 
actly identical to the energy gain of the system in our 
simulations. 

3.2.4. Energy Dissipation at the Nonlinear Stage 

The role of dissipation processes in the energy con- 
version is investigated in this subsection. The injected 
energy (Ein) — {Eh-i.w) + {E,n,v) should be balanced with 
the increase of the sum of Eth — P/{j—l), E^n — B'^/8n, 
and E'k = pv'^/2. Figure [6^ depicts the time evolution 
of (ii'th), {Em), and (E^) for the case with i?MRi — 0.01. 
The viscous heating rate (Evis) is also shown in this fig- 
ure, which is defined as 



E^ 



(17) 
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where $ij is the same definition that is used in the energy 
equation (fTTj) . 

The sum of the components perfectly coincides with 
the energy gain of the system. Using the time- and 
volume-averaged values, which are indicated by double 
brackets, the gain rates of each energy component are 
{{Ea,))/m^)) = 0.233, ((E,^)) / ((Ei^)) = 0.757, and 
{{E-k})/ {{Eiu)) ^ 0.01. Here the time average is taken at 
18 < i/^rot < 25. This indicates that a large portion of 
the injected energy is converted into the magnetic energy 
of the system. 

Because the system is almost incompressible, the dom- 
inant heating mechanism in our simulations should be 
the kinematic viscous heating. However the thermal en- 
ergy gain is much larger than the viscous heating rate 
(see Fig.|B^). The ratio of these components is evaluated 

as ((£;vis))/((£^th)) = 0.36 in the range 18 < t/trot < 25. 
This fraction becomes smaller and smaller for the cases 
with larger Reynolds number. The rest of the heating is 
caused by the numerical magnetic dissipation. Then our 
results suggest that the magnetic dissipation is the dom- 
inant mechanism of heating at the nonlinear phase of the 
MRI and might play an essential role for the suppression 
of MRI driven turbulence. 

Actually in the resistive MHD cases, the heating of the 
system is controlled by the joule heating almost com- 
pletely (Sano & Inutsuka 2001). The energy changing 
rates in a very resistive model with the Lundquist num- 
ber ^MRi = 0.1 is demonstrated in Figure [61d. This is a 
typical case in which the channel flow is disrupted and 
MHD turbulence is damped at the nonlinear stage (see 
Fig. [4^). The ratio of the total energy gain and the ther- 
mal energy gain is {{Eth)) / {{Ein)) — 0.99 in the range 
20 < t/trot < 25. We find that the joule heating takes a 
major role in the thermal energy gain of the system, that 

is {{E^ou)) / {{Eth)) — 0.93. In the resistive system, the 
magnetic energy amplified by the MRI is transformed 
into the thermal energy via joule heating or magnetic 
reconnection throughout the nonlinear evolution. 

3.2.5. Stress at the Nonlinear Stage 

Finally, the time- and volume-averaged atot at the 
nonlinear stage is depicted as a function of the initial 
Reynolds number i?MRi and Lundquist number S'mri in 
Figure[71 We take the average of (atot) over 5 orbits just 
after the time when the ratio (£'th) / {Em) begins to rise 
and the nonlinear evolution is started. Diamond-shape 
shows the results in the viscous fiuid and cross-shape is 
that in the resistive one. Note that these are not the sat- 
urated values. Upward arrow over-plotted on the sym- 
bols denotes that the value is the lower limit, because 
((o^tot)) is still increasing with time. The downward ar- 
row stands for decaying models and thus the stress is the 
upper limit. 

The stress at the nonlinear stage are almost the same 
when the diffusion is weak (-Rmri ^ 1 or 5mri ^ 1)- 
However, a huge difference can be seen in the highly dif- 
fusive regime. In the presence of the ohmic dissipation, 
the stress rapidly decreases with decreasing S'mri- For 
the models with the kinematic viscosity, on the other 
hand, it increases with the decrease of i?MRi- The origin 
of this difference between viscous and ohmic dissipative 
systems is discussed later in § 4.1. 



The inverse correlation between ((atot)) and -Rmri for 
the cases with large viscosity could be originated from 
stable growth of a channel flow. This is because the 
large viscosity can suppress the growth of any other 
modes than the two-channel flow, and thus the chan- 
nel mode can evolve up to highly nonlinear amplitude. 
The viscosity may enhance the saturation amplitude of 
the MRI. However, our results are restricted in two- 
dimensional simulations. The nonlinear evolution of the 
MRI in three-dimension must be quite different, because 
the channel ffow is known to be unstable to the nonax- 
isymmetric parasitic instability (Goodman & Xu 1994). 
We are planning to perform three-dimensional study of 
the MRI for verifying these nonlinear properties in the 
viscous accretion disks. 

4. DISCUSSION 
4.1. Nonlinear Behavior in the Single Diffusive Systems 

To give a physical explanation for the nonlinear be- 
havior of the axisymmetric MRI, we focus on the critical 
wavelength obtained from the linear theory in this sec- 
tion. The diagrams of Figures [Ha, and [Sh indicate the 
critical and the fastest growing wavelengths of the MRI 
as a function of the Lundquist number Smki and the 
Reynolds number i?MRi, respectively (see also Table [l}. 
Note that the vertical axes are normalized by 2TT{r]/V.)^^^ 
in Figure Hi and by 2n{v/n)^^^ in Figure [HId. Shaded 
area denotes the linearly unstable regions for the MRI. 
Assuming fixed diffusivities, S'mri and i?MRi increase as 
the instability grows because they are proportional to 
the squared Alfven velocity. Then the horizontal axis in 
Figure [8] can be regarded as the time direction in terms 
of the evolution of the MRI. 

First, we consider the resistive case shown by Fig- 
ure [8^. For the case of Smri £ 1, the critical wavelength 
is described as Adit — 'hIva- At the nonlinear stage of 
the two-dimensional MRI, MHD turbulence decays and 
it saturates only when Smri ^ 1 (see Figs. SJd and[71). 
This behavior can be interpreted schematically using the 
Acrit-SMRi diagram (Sano & Miyama 1999). As the MRI 
grows and amplifies the magnetic field, the critical wave- 
length shifts to the shorter length-scale. Then, many 
smaller scale fluctuations can become unstable. Those 
structures enhance the efficiency of ohmic dissipation in 
the turbulent state. In other words, the system evolves 
toward a more dissipative state, and could be saturated 
at a critical point around Smri — 1, at which the criti- 
cal wavelength reaches the shortest value and the ohmic 
dissipation is the most efficient. In this way, MHD tur- 
bulence can decay if Smri ^ 1- 

When Smri ?i 1; on the other hand, the critical wave- 
length is given by Acrit — vaI^- Two-dimensional cal- 
culations of the MRI suggests that the unstable growth 
cannot saturate if Smri ^ 1- The nonlinear behavior in 
these cases can be explained by Figure [5^ analogously 
as follows: At the linear evolutionary stage, the critical 
wavelength in the radial direction becomes longer due 
to the exponential growth of the radial field component, 
while the vertical field grows slower than exponentially. 
The wavevectors of the unstable modes thus become par- 
allel to the vertical axis. This channel flow mode is the 
exact solution of the nonlinear MHD equations (Good- 
man & Xu 1994). As the field strength becomes larger. 
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the critical wavelength shifts to the longer one. The dis- 
sipation can be much less effective, and thus the channel 
solution continues to grow without saturation. 

Next, let us consider the viscous case shown in Fig- 
ure [HId. In the viscous fluid, the critical wavelength is 
given by Acrit — va/^ despite the size of the Reynolds 
number (see Figure [T]) . Even if i?MRi is much smaller 
than unity, the critical wavelength thus shifts to larger 
scale as the instability grows. Then the system always 
evolves toward a less dissipative state and is not satu- 
rated. This interpretation is consistent with our numer- 
ical results shown in Figure [T] 

These results indicates that the saturation process of 
the MRI would be changed dramatically at the critical 
point at which the critical wavelength switches from the 
decreasing function of the field strength to the increasing 
one. The differences in the nonlinear behavior of the 
MRI between the viscous and resistive systems would be 
originated from whether the critical point exists or not. 

4.2. MRI in the Doubly Diffusive System 

In this subsection, we apply the discussion above to the 
doubly diffusive system that includes both the viscosity 
and resistivity. The saturation behavior of the axisym- 
metric MRI can be anticipated by the dependence of the 
critical wavelength on the field strength derived from the 
linear theory. With the same framework used in § 2, a 
local axisymmetric dispersion equation of the MRI in the 
presence of both viscous and ohmic dissipations is given 

by 

047'* -I- 037^ -I- a27^ -I- ai7 oq = , (18) 

where 

04 = 1 , 03 = 2 ( h ) kl , 

V ^MRI JMKl J 

02 = ( + ^ + -f^^ ) ~kl + 2kl + , 

(-Rmri •S'mri) -Rmri^mri 

V -KmRI OMRI / OMRI 

nn - ^ P 4- ^ /~.6 

°0 ~ n2 c2 + z? c '^z 

^mri'^mri ^mri'-'mri 

\ '-'mri / 

(Menou et al. 2006; Masada et al. 2007; Lesur & Lon- 
garetti 2007; Pessah & Chan 2008). This equation is, as 
expected, characterized by i?MRi and S'mri- 

Here we focus on the system with a constant mag- 
netic Prandtl number (Pm = Sum/ Rum — v/rf). Fig- 
ure [9] demonstrates the critical wavelength of the MRI as 
a function of S'mri for various values of Pm obtained 
by solving the dispersion equation (18). The 5mri- 
dependence of the critical wavelength varies with the 
size of Pm. When Pm ^ 1, the linear growth of the 
MRI is independent of the magnetic Prandtl number, 
and the critical wavelength is almost identical to the pure 
resistive case [Pm = 0). However, if the viscosity ef- 
fects is added sufficiently, then the critical wavelength 
is enlarged by the suppression due to the viscosity in 



the middle range of S'mri. For the cases of Pm ^ 1, 
the critical wavelength around Smri — 1 is given by 
k-^l cx (SmriPmri)-^/^ « {Sl,^^/Pm)-^/^ (Pessah & 
Chan 2008). 

Although the critical wavelength shifts to longer 
length-scales as the magnetic Prandtl number increases, 
the critical wavelength has a minimum value for all the 
cases. This implies that the MRI turbulence could be 
suppressed if the Lundquist number is less than a critical 
value. This diagram suggests that the critical Lundquist 
number S'mri, c depends on Pm. The critical Lundquist 
number S'mri, c is plotted as a function of Pm in Fig- 
ure [10^. In the regime of Pm ^ 1, it is proportional to 
the square root of the magnetic Prandtl number, that is 
Scrit,c oc Pm^/'^. In contrast, it remains to be constant 
in the range Pm <C 1. 

Since the critical Reynolds number is given by 
Rmri,c — Smki,c/ Pfn, we can also obtain the relation 
between Pm and Pmri,c, and which is depicted in Fig- 
urefTOb. Nonlinear growth of the MRI can be expected in 
the parameter region above this critical curve. The mag- 
netic Prandtl number is proportional to Pm|^j ^ in the 

regime of Pm 3> 1 and Pm cx Pm^j ^ when Pm ^ 1 . 
This curve is reminiscent of the critical curve for MHD 
turbulence sketched from nonlinear simulations of MRI 
(Fromang et al. 2007). The Reynolds number Pmri 
in their models are at most a few tens^, and the critical 
magnetic Prandtl number is around unity. Thus our pre- 
diction is roughly consistent with nonlinear results even 
quantitatively. Note that the critical curve shown by 
Figure [10] is obtained by using only the features of the 
linear dispersion relation of MRI. This implies that the 
linear growth of the MRI is very important even at the 
nonlinear saturated phase to sustain MHD turbulence. 

The discussion in this paper is based only on two- 
dimensional simulations of the MRI. For the understand- 
ing the saturation mechanism of the MRI, it is quite 
important to perform the systematic three-dimensional 
analysis of the MRI in the presence of multiple diffusiv- 
ities. Furthermore, the assumption of the local shearing 
box could affect the nonlinear evolution of the viscous 
MRI. The necessary ingredients of the unstable growth of 
the MRI are the velocity shear and the magnetic field. In 
the numerical setting of the local shearing box, the veloc- 
ity profile of the background shear flow cannot disappear 
by the role of the kinematic viscosity, but is imposed by 
the boundary conditions. It would be very interesting to 
use global disk models to investigate the MRI in highly 
viscous disks. These are our next tasks. 

5. SUMMARY 

Axisymmetric MRI in viscous accretion disks is investi- 
gated by linear and nonlinear analyses. A local shearing 
box threaded by a uniform vertical magnetic field is used 
for our nonlinear simulations. The nonlinear results of 
the viscous MRI are compared with the resistive case fo- 
cusing on two non-dimensional parameters, the Reynolds 
number for the MRI Pmri and the Lundquist number 
for the MRI Smri- Our main findings are summarized 

The definition of the Reynolds number Re in Fromang et al. 
(2007) is different from -Rmri in this paper. The relation between 
these two is -Rmri ~ ctM^e, where ojv/ is the Maxwell stress nor- 
malized by the (initial) pressure. 
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as follows. 

1. In axisymmetric two-dimensional simulations, the 
MRI continues growing regardless of the size of the 
Reynolds number. When i?MRi ^ 1, the stress in its 
nonlinear stage is inversely correlated with i?MRi, and 
thus can be larger than that in the ideal MHD run. In 
the highly resistive fluid, on the other hand, the growth 
of the MRI is saturated and MHD turbulence dies away. 
When the Lundquist number is less than unity, the satu- 
rated stress decreases dramatically with decreasing S'mri- 

2. At the nonlinear stage of the MRI, a large portion 
of the injected energy is converted to the magnetic en- 
ergy in the viscous system (((£;,„)) > {{Eth}) > ((i^k)))- 
The thermal energy gain is much larger than the viscous 
heating rate. In contrast, for the case of the resistive 
system, the thermal energy is converted from the mag- 
netic energy through the joule heating. This difference 
in the energy dissipation efficiency may affect the satu- 
ration process of the MRI. 

3. Nonlinear behavior of the MRI in the single diffu- 



sive system can be understood with the help of the local 
dispersion relation. The key characteristic is the depen- 
dence of the critical wavelength on the Reynolds number 
-Rmri and the Lundquist number 5'mri- Applying this 
interpretation to the doubly diffusive system with both 
the viscous and ohmic dissipations, a condition for sus- 
taining MRI driven turbulence is obtained as a function 
of i?MRi and S'mri- 
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Ideal MHD 


Resistive Case 


Viscous Case 


Balancing Rate 


Alfvcn frequency 


Alfvcn frequency 


Shear rate 




= Rotation frequency 


= Dissipation rate 


= Dissipation rate 










Unstable Wavelength [A = fc~i] 


X = va/^ 


A = ri/vA 


A = (i//a)i/2 


Growth Rate [7 = va/A] 


7 = n 


7 = v\/r) 


7 = (t;in/i/)V2 



TABLE 1 

Unstable wavelength and growth rate op the MRI for ideal MHD, resistive, and viscous cases. 
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Fig. 1. — Linear growth rate of the MRI as a function of the wavenumbcr. The cases with different Reynolds numbers for MRI 
-Rmri = v^/uQ = oo, 10, 1.0 and 0.1 are depicted. Normalization of the vertical and horizontal axes are the angular velocity H and the 
typical wavenumbcr of the MRI va/^, respectively. The symbols plotted over the analytical dispersion relation are numerical growth rates 
calculated by our numerical scheme. 
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Fig. 3. — Time evolution of the volume-averaged (a) Reynolds stress {ocr) = {wii)/{P) and (b) Maxwell stress (om) = {wm)/{P) for 
the cases with different Reynolds number Rmri = 10.0, 1.0, 0.1 and 0.01. The horizontal axis is normalized by the disk rotation time 
trot = 27r/f2. 
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Fig. 4. — Panel (a): Time evolution of the volume-averaged a parameter of Shakura & Sunyaev (1973) in the viscous fluid for the cases 
with different Reynolds number -Rmri = 10.0, 1.0, 0.1 and 0.01. Panel (b): The a parameter in the resistive fluid for the cases with 
different Lundquist number Smri = 100.0, 10.0, 1.0, 0.3 and 0.1. Normalization of the horizontal axis is the same as Figure|3] 
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Fig. 5. — Time evolution of the volumc-avcragcd time derivative of the total energy (Etot) = (-S'th) + {Em) + (Eu) and the input 

energies due to the turbulent stress {-Ein.-m} and the kinematic viscous stress (-Ein,t))- These quantities should satisfy the energy conservation 

(^'tot) = {Eir\,w) + {Eiii,v) ■ This is for the case with i?MRi = 0.01. The vertical axis is given in the unit of Etho/^rot where i?thO = -Po/(7 — 1) 
is the initial thermal energy. 
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Fig. 6. — Time evolution of the volume- averaged time derivative of the input energy (Ein) = (^in,™) + (Ein^v), thermal energy {Eth)y 
magnetic energy (Em), kinetic energy (Ek)- Panel (a) shows the result of the viscous case with -Rmri = 0.01 in the period 16 < t/trot < 25. 
The volume- averaged viscous heating rate (Eyis) is shown in this figure. Panel (b) depicts that of the resistive case with Smri = 0.1 in the 



period 20 < t/Uot < 25. 
as those in Figure |4] 



The volume- averaged joule heating rate (Ejou) is shown in this figure. Normalizations of each axis are the same 
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Fig. 7. — Time- and volume-averaged atot at the nonlinear stage as a function of the initial Reynolds number i?MRl and Lundquist 
number Smri- We take the time-average of (atot) over 5 orbits at the nonlinear regime. Diamonds are the results in the viscous fluid and 
crosses are those in the resistive one. Note that these are not the saturated values. The upward arrow denotes models in which ((atot)) is 
still increasing with time and the downward arrow stands for decaying models. 
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Fig. 8. — Characteristic wavelengths of the MRI for (a) the resistive case and (b) the viscous case. The horizontal axis in each panel is 
(a) the Lundquist number 5mri and (b) the Reynolds number i?MRl- Shaded area denotes the unstable regions for the MRI expected from 
the linear theory. The critical wavelength in the resistive case takes a minimum value, while it monotonically increases in the viscous case. 
The diflference in the nonlineaj: regime between the resistive and viscous models can be explained by this feature in the critical wavelength 
of the MRI. 
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Fig. 9. — Schematic picture of the critical wavelength of the MRI as a function of the Lundquist number S'mri for different values of 
the magnetic Prandtl number. Thick black curve represents the critical wavelength for the cases with Pm <S 1. The models with Pm ^ 1 
are depicted by yellow (Pm = 10), blue (Pm = 10^), green (Pm = 10^), and orange (Pm = 10^) curves. Since the diffusive parameters 
are fixed, the Lundquist number S'mri is a function of the field strength. Shaded area denotes the unstable regions for the MRI expected 
from the linear theory at small Prandtl numbers. The critical point and critical Lundquist number are marked by the filled red and black 
circles, respectively. 
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Fig. 10. — (a) The critical Lundquist number as a function of the magnetic Prandtl number Pm. (b) T he magnetic Prandtl number as a 
function of the critical Reynolds number. These relations are derived by solving the dispersion equation I II8I I. The parameter region above 
the critical curve denotes where the nonlinear growth of the MRI can be expected and MRI driven turbulence will be sustained. 



